home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / regress.pro < prev    next >
Text File  |  1997-07-08  |  5KB  |  130 lines
  1. ; $Id: regress.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1982-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. FUNCTION REGRESS,X,Y,Weights,YFIT,Const,SIGMA,FTEST,R,RMUL,CHISQ,STATUS, $
  7.                  RELATIVE_weight=relative_weight
  8. ;+
  9. ; NAME:
  10. ;    REGRESS
  11. ;
  12. ; PURPOSE:
  13. ;    Perform a multiple linear regression fit.
  14. ;
  15. ;    REGRESS fits the function:
  16. ;        Y[i] = Const + A[0]*X[0,i] + A[1]*X[1,i] + ... + 
  17. ;                      A[Nterms-1]*X[Nterms-1,i]
  18. ;
  19. ; CATEGORY:
  20. ;       G2 - Correlation and regression analysis.
  21. ;
  22. ; CALLING SEQUENCE:
  23. ;    Result = REGRESS(X, Y, Weights, Yfit, Const, Sigma, Ftest, R, Rmul, Chisq)
  24. ;
  25. ; INPUTS:
  26. ;       X:    The array of independent variable data.  X must 
  27. ;        be dimensioned as an array of Nterms by Npoints, where 
  28. ;        there are Nterms coefficients (independent variables) to be 
  29. ;        found and Npoints of samples.
  30. ;
  31. ;       Y:    The vector of dependent variable points.  Y must have Npoints 
  32. ;        elements.
  33. ;
  34. ; Weights:    The vector of weights for each equation.  Weights must be a vector
  35. ;        of Npoints elements.  For instrumental (Gaussian) weighting, 
  36. ;        W[i] = 1/standard_deviation(Y[i])^2.  For statistical  (Poisson)
  37. ;        weighting, w[i] = 1./Y[i].  For no weighting, set w[i]=1,
  38. ;        and also set the RELATIVE_WEIGHT keyword.
  39. ;
  40. ; OUTPUTS:
  41. ;    REGRESS returns a column vector of coefficients that has Nterms 
  42. ;    elements.
  43. ;
  44. ; OPTIONAL OUTPUT PARAMETERS:
  45. ;    Yfit:    Vector of calculated values of Y with Npoints elements.
  46. ;
  47. ;      Const:    Constant term. (A0)
  48. ;
  49. ;    Sigma:    Vector of standard deviations for coefficients.
  50. ;
  51. ;    Ftest:    The value of F for test of fit.
  52. ;
  53. ;    Rmul:    The multiple linear correlation coefficient.
  54. ;
  55. ;    R:    Vector of linear correlation coefficients.
  56. ;
  57. ;    Rmul:   The multiple linear correlation coefficient.
  58. ;
  59. ;    Chisq:    Reduced, weighted chi squared.
  60. ;
  61. ;       Status:  A named variable to receive the status of the INVERT 
  62. ;                (array inversion) computation. A value of 0 indicates 
  63. ;                a successful computation. A value of 1 indicates the 
  64. ;                inversion is invalid due to a singular array. A value 
  65. ;                of 2 indicates the possibility of an inaccurate result 
  66. ;                due to the use of a small pivot element.
  67. ;
  68. ; KEYWORDS:
  69. ; RELATIVE_WEIGHT:  If this keyword is set, the input weights
  70. ;        (W vector) are assumed to be relative values, and not based
  71. ;        on known uncertainties in the Y vector.  Set this keyword in 
  72. ;        the case of no weighting, W[*] = 1.
  73. ;
  74. ; PROCEDURE:
  75. ;    Adapted from the program REGRES, Page 172, 
  76. ;    Bevington, Data Reduction and Error Analysis for the 
  77. ;    Physical Sciences, 1969.
  78. ;
  79. ; MODIFICATION HISTORY:
  80. ;    Written, DMS, RSI, September, 1982.
  81. ;    Added RELATIVE_WEIGHT keyword    W. Landsman   August 1991
  82. ;       Fixed bug in invert  Bobby Candey 1991 April 22
  83. ;       Added STATUS argument.  GGS, RSI, August 1996
  84. ;-
  85. ;
  86. On_error,2              ;Return to caller if an error occurs 
  87. SY = SIZE(Y)            ;Get dimensions of x and y.  
  88. SX = SIZE(X)
  89. IF (N_ELEMENTS(Weights) NE SY[1]) OR (SX[0] NE 2) OR (SY[1] NE SX[2]) THEN $
  90.   message, 'Incompatible arrays.'
  91. ;
  92. NTERM = SX[1]           ;# OF TERMS
  93. NPTS = SY[1]            ;# OF OBSERVATIONS
  94. ;
  95. SW = TOTAL(Weights)           ;SUM OF WEIGHTS
  96. YMEAN = TOTAL(Y*Weights)/SW   ;Y MEAN
  97. XMEAN = (X * (REPLICATE(1.,NTERM) # Weights)) # REPLICATE(1./SW,NPTS)
  98. WMEAN = SW/NPTS
  99. WW = Weights/WMEAN
  100. ;
  101. NFREE = NPTS-1          ;DEGS OF FREEDOM
  102. SIGMAY = SQRT(TOTAL(WW * (Y-YMEAN)^2)/NFREE) ;Weights*(Y(I)-YMEAN)
  103. XX = X- XMEAN # REPLICATE(1.,NPTS)      ;X(J,I) - XMEAN(I)
  104. WX = REPLICATE(1.,NTERM) # WW * XX      ;Weights(I)*(X(J,I)-XMEAN(I))
  105. SIGMAX = SQRT( XX*WX # REPLICATE(1./NFREE,NPTS)) ;Weights(I)*(X(J,I)-XM)*(X(K,I)-XM)
  106. R = WX #(Y - YMEAN) / (SIGMAX * SIGMAY * NFREE)
  107. ARRAY = (WX # TRANSPOSE(XX))/(NFREE * SIGMAX #SIGMAX)
  108. IF (SX[1] EQ 1) THEN ARRAY = 1 / ARRAY ELSE begin
  109.   ARRAY = INVERT(Array, status)
  110.   if status eq 1L then MESSAGE, "Inversion Failed due to singular array."  
  111.   endelse
  112. A = (R # ARRAY)*(SIGMAY/SIGMAX)         ;GET COEFFICIENTS
  113. YFIT = A # X                            ;COMPUTE FIT
  114. Const = YMEAN - TOTAL(A*XMEAN)             ;CONSTANT TERM
  115. YFIT = YFIT + Const                        ;ADD IT IN
  116. FREEN = NPTS-NTERM-1 > 1                ;DEGS OF FREEDOM, AT LEAST 1.
  117. CHISQ = TOTAL(WW*(Y-YFIT)^2)*WMEAN/FREEN ;WEIGHTED CHI SQUARED
  118. ;
  119. ; If all the weights are 1 then
  120. ; chisq = variance 
  121. ;
  122. IF KEYWORD_SET(relative_weight) then varnce = chisq $
  123.                                 else varnce = 1./wmean
  124. sigma = sqrt(array[indgen(nterm)*(nterm+1)]*varnce/(nfree*sigmax^2)) ;Error term
  125. RMUL = TOTAL(A*R*SIGMAX/SIGMAY)         ;MULTIPLE LIN REG COEFF
  126. IF RMUL LT 1. THEN FTEST = RMUL/NTERM / ((1.-RMUL)/FREEN) ELSE FTEST=1.E6
  127. RMUL = SQRT(RMUL)
  128. RETURN,A
  129. END
  130.